###此代码用于分析727个位点能否富集到某个组织，组织特异性的result用的是师兄算好的文件
library(meta)
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/")
file=read.csv("20201120/at.least.one.AShM.in.DC.add.BF.beta0.add.CCHC.csv",head=T)
filea=read.csv("20201112做汇总表/all.FDR.sig.at.least.one.add.direction.same.diff.csv",head=T)
filea$id=paste(filea$Chr,filea$Start,sep = ":")
filea1=filea[filea$FDR.sig>1,]

file$id=paste(file$Chr,file$Start,sep=":")
file1=file[file$pattern.not.rm.dupl.num.DC>1,]
file2=file1[file1$BF_in_DC>1,]
file3=file1[file1$BF_in_DC>10,]


tissue.spec=list.files(pattern = "txt",path="E:/5hmc_file/组织特异性表达/",recursive = T,full.names = T)##加入了递归显示，方便提取子目录下的文件

#case为807，背景为随机150K时，几乎所有p值显著的OR<1，因此这个组合不行
#case为807，背景为117K时，p值显著的OR<1，因此这个组合不行
#case为117K，背景为随机150K时,p值显著的OR<1，因此这个组合不行

con1=read.csv("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201207/con_genotype.random1.anno.hg19_multianno.csv",head=T)
con1.gene=unique(unlist(strsplit(as.character(con1$Gene.refGene),";")))

con.file=list.files(pattern="con_genotype.random",path="./20201207/",full.names=T)

case=filea
caseid=unique(unlist(strsplit(as.character(case$Gene.refGene),";")))
caseid=caseid[!caseid=="NONE"]

col_names=c("tissue","spec_filename","con.file","case_overlap","case_not","con_overlap","con_not","OR","p.value")
result=data.frame(matrix(NA,1,ncol=9))
names(result)=col_names
for(l in 1:length(con.file)){
con=read.csv(con.file[l],header=T)
con.gene=unique(unlist(strsplit(as.character(con$Gene.refGene),";")))
con.gene=con.gene[!con.gene=="NONE"]

for(k in 1:length(tissue.spec)){
file=read.table(tissue.spec[k],head=T,sep="\t")###筛选出有组织特异性的基因
file$num=rowSums(file[,3:ncol(file)]>0.9,na.rm=T)
file1=file[file$num>0,]
cn=names(file1)


for (j in 3:(ncol(file1)-1)) {
  tmp=data.frame(Description=file1$Description,file1[,j],file1[,j])
  tmp$num=rowSums(tmp[,2:3]>0.9,na.rm = T)
  sel=which(tmp$num>1)
  tmpid=tmp[sel,]$Description
  result_tmp=data.frame(matrix(NA,1,ncol=9))
  names(result_tmp)=col_names
  result_tmp[,1]=gsub("_spec","",cn[j])
  result_tmp[,2]=gsub("E:/5hmc_fil/组织特异性表达/","",gsub("./组织特异性表达/","",tissue.spec[k]))
  result_tmp[,3]=gsub("./20201207/","",gsub(".hg19_multianno.csv","",con.file[l]))
  result_tmp[,4]=length(intersect(caseid,tmpid))
  result_tmp[,5]=length(caseid)-length(intersect(caseid,tmpid))
  result_tmp[,6]=length(intersect(con.gene,tmpid))
  result_tmp[,7]=length(con.gene)-length(intersect(con.gene,tmpid))
  result_tmp[,8]=fisher.test(matrix(c(result_tmp[1,4],result_tmp[1,5],result_tmp[1,6],result_tmp[1,7]),nrow = 2))$estimate
  result_tmp[,9]=fisher.test(matrix(c(result_tmp[1,4],result_tmp[1,5],result_tmp[1,6],result_tmp[1,7]),nrow = 2))$p.value
  result=rbind(result,result_tmp)
}
}
}
result=result[-1,]
result1=result[result$p.value<0.05,]
write.csv(result,"tissue_spec_enrichment_result.csv",quote=F,row.names = F)
write.csv(result1,"tissue_spec_enrichment_result_sig.csv",quote=F,row.names = F)

result_c2=result1[result1$con.file=="nobias_AShM",]
result_c2=result_c2[result_c2$OR>1,]
result_c2=result_c2[order(result_c2$OR),]
result_c2$case_not=result_c2$case_overlap+result_c2$case_not
result_c2$con_not=result_c2$con_overlap+result_c2$con_not

metaor3<-metabin(case_overlap,case_not,con_overlap,con_not,data=result_c2,sm="OR",studlab = tissue) 
forest(metaor3)

dup_id=unique(result_c2[duplicated(result_c2$tissue),]$tissue)###把重复行挑出来重新处理求均值,然后求OR和pvalue
result_c2_no_dup=result_c2[!result_c2$tissue %in% dup_id,]
for(m in 1:length(dup_id)){
rt_tmp =result_c2[result_c2$tissue==dup_id[m],]
rt_tmp2=rt_tmp[1,]
rt_tmp2$spec_filename=paste(unlist(gsub(".txt","",unique(rt_tmp$spec_filename))),collapse = "_")
rt_tmp2$case_overlap=round(mean(rt_tmp$case_overlap))
rt_tmp2$con_overlap=round(mean(rt_tmp$con_overlap))
rt_tmp2[,8]=fisher.test(matrix(c(rt_tmp2[1,4],rt_tmp2[1,5],rt_tmp2[1,6],rt_tmp2[1,7]),nrow = 2))$estimate
rt_tmp2[,9]=fisher.test(matrix(c(rt_tmp2[1,4],rt_tmp2[1,5],rt_tmp2[1,6],rt_tmp2[1,7]),nrow = 2))$p.value
result_c2_no_dup=rbind(result_c2_no_dup,rt_tmp2)
}

result_c2_no_dup=result_c2_no_dup[order(result_c2_no_dup$OR),]
metaor_no_dup<-metabin(case_overlap,case_not,con_overlap,con_not,data=result_c2_no_dup,sm="OR",studlab = tissue) 
forest(metaor_no_dup)										######森林图
															######erro bar的图
meta_data=as.data.frame(metaor_no_dup)
meta_data$issue=gsub("Brain...","",result_c2_no_dup$tissue)
meta_data=meta_data[order(meta_data$TE),]
meta_data$num=1:dim(meta_data)[1]
ggplot(data=meta_data,aes(y=TE,x=num,color=pval))+geom_point(size=3)+
    geom_errorbar(aes(ymin=lower,ymax=upper),width=0.1,color="black")+coord_flip()+scale_color_gradient(low = "red", high = "green")+
    scale_x_continuous(breaks = meta_data$num,labels = meta_data$issue)+theme_light()+ylab("log(OR)")+
    theme(panel.grid.minor = element_blank())
	
															######有标记的散点图
library(ggrepel)
meta_data1=data.frame(issue=meta_data$issue,TE=meta_data$TE,num=meta_data$num)
 meta_data2=meta_data1
 meta_data2$TE=0
 meta_data_line=rbind(meta_data1,meta_data2)
ggplot(meta_data,aes(TE,num))+geom_point(aes(size=event.e,color=pval))+geom_text_repel(aes(label = event.e))+
  scale_color_gradient(low = "red", high = "green")+geom_line(data = meta_data_line,aes(x=TE,y=num,group=issue))+
  labs(x="log(OR)",y="issue")+scale_y_continuous(breaks = meta_data$num,labels = meta_data$issue)+
  theme_bw()+theme(panel.grid.minor = element_blank(),panel.grid.major = element_blank())